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Abstract 

Evolutions of the trading landscape lead to the capability to exchange the same financial 
instrument on different venues. Because of liquidity issues, the trading firms split large orders 
across several trading destinations to optimize their execution. To solve this problem we devised 
two stochastic recursive learning procedures which adjust the proportions of the order to be sent 
to the different venues, one based on an optimization principle, the other on some reinforcement 
ideas. Both procedures are investigated from a theoretical point of view: we prove a.s. conver- 
gence of the optimization algorithm under some light ergodic (or "averaging") assumption on 
the input data process. No Markov property is needed. When the inputs are i.i.d. we show that 
the convergence rate is ruled by a Central Limit Theorem. Finally, the mutual performances of 
both algorithms are compared on simulated and real data with respect to an "oracle" strategy 
devised by an "insider" who a priori knows the executed quantities by every venues. 

Keywords Asset allocation, Stochastic Lagrangian algorithm, reinforcement principle, monotone 
dynamic system. 

2001 AMS classification: 62L20, secondary: 91B32, 62P05 

1 Introduction 

The trading landscape have seen a large number of evolutions following two regulations: Reg 
NMS in the US and MiFID in Europe. One of their consequences is the capability to exchange 
the same financial instrument on different trading venues. New trading destinations appeared to 
complement the trading capability of primary markets as the NASDAQ and the NYSE in the US, 
or EURONEXT, the London Stock Exchange and Xetra in Europe. Such alternative venues are 
called "Electronic Communication Network" (ECN) in the US and Multilateral Trading Facilities 
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(MTF) in Europe. Each trading venue differentiates from the others at any time because of the 
fees or rebate it demands to trade and the liquidity it offers. 

As the concerns about consuming liquidity increased with the financial crisis, trading firms use 
Smart Order Routers (SOR) as a key element in the process of optimizing their execution of large 
orders. Such devices are dedicated to split orders across trading destinations as a complement to 
the temporal slicing coming from the well known balance between the need to trade rapidly (to 
minimize market risk) and trading slow (to avoid market impact). 

If the temporal slicing has been studied since the end of the nineties pQ with recent advances 
to adapt it to sophisticated investment strategies [19], this kind of spatial slicing (across trading 
destinations) has been mainly studied by economists from the point of view of its global market 
efficiency [8] rather than from one investor's point of view. 

The complexity of spreading an order between N trading destinations comes from the fact that 
you never knows the quantity Di available on the i th trading venue to execute your order of size 
V during a time interval dt at your given price. If the fraction r{ V of your order that you sent to 
the i th liquidity pool is higher than D^. you will loose time and may loose opportunity to execute 
T{ V — D{ in an another pool; on another hand if r, V is lower than D^. you will loose money if this 
pool fees are cheap, and an opportunity to execute more quantity here. The only way to optimize 
such a split on real time is to adjust on the fly the proportions (rj)j according to the result of your 
previous executions. 

This paper is an in depth analysis of the optimal split of orders. The illustrations and most 
of the vocabulary come from the "Dark pool" case, where the price S is not chosen by the trader 
(it is the market "mid point" price) and the answer of the pool is immediate (i.e. St = 0). Dark 
pools are MTFs that do not publish pre-trade informations, so an efficient use of the results of 
the previous executions (namely the realizations of the min(Z?|, r\ V ) for any i and all t in the 
past) is crucial. The results exposed here solve the problem of simultaneously splitting orders and 
using the information coming back from the pools to adjust the proportions to send for the next 
order, according to a criteria linked to the overall quantity executed (i.e. a linear combination of 
the min(Di,riV)). 

The resulting trading strategy (which optimality is proven here) can be considered as an exten- 
sion of the one conjectured by Almgren in [2]. It may also be related to the class of multi-armed 
bandit recursive learning procedures, recently brought back to light in several papers (see [HI [22], 
|15l [T6] ; which in turn belongs to the wide family of "recursive stochastic algorithms" also known 
as "stochastic approximation" and extensively investigated in the applied probability literature 
(see m, M, U\, etc)). 

In fact, we introduce two learning algorithms one based on an optimization under constraints 
principle and a second algorithm based on a reinforcement principle for which we establish the 
existence of an equilibrium. We extensively investigate the first one, considering successively the 
classical - although unrealistic - case where the inputs (requests, answers) are i.i.d. and a setting 
in which the input only share some averaging properties. In the i.i.d. setting we establish a.s. con- 
vergence of the procedure and a Central Limit Theorem relying on classical results from Stochastic 
Approximation Theory. By averaging setting (also referred as ergodic setting), we mean that the 
inputs of the procedure has a.s. an averaging property with respect to a distribution v at a given 
rate, say n~@, f3 > 0, for a wide enough class of Borel functions. Typically, in our problem, these 
inputs are the successive ./V + 1-tuples (V n , Df, i = 1, . . . , N), n > 1. Typically, if we denote this 
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input sequence inputs by (Y n ) n >i, we will assume that, for every f £ V/3 ;P , 

n ,. 

- V f(Y k ) - fdv = Ofn"' 3 ) P-a.s and in L P (P). 

n ^ 1 

Usually, V^ iP is supposed to contain at least bounded continuous function g : R^ +1 — > R and 
subsequently all bounded u-a.s. continuous functions. This will be enough for our purpose in this 
paper (Stochastic approximation in this general framework is investigate in |17j). But the key point 
to be noted here is that no Markov assumption is needed on this input sequence (Y n ) n >i. These 
assumptions are hopefully light enough to be satisfied by real data since it can be seen as a kind 
of "light" ergodicity at a given rate. In a Markovian framework it could be related to the notion 
of "stability" in the literature, see [Tj- 

Thus, this setting includes stationary a-mixing processes (satisfying an Ibragimov condition) 
like those investigated in [6] (in [5] weaker dependence assumptions are made in the chapter devoted 
to stochastic approximation but the perturbation is supposed to be additive and non causal which 
is not at all the case in our problem). As concerns the second procedure for which no Lyapunov 
function seems to be (easily) made available, we establish the existence of an equilibrium and 
show the ODE related to the algorithm is a competitive system in the terminology of monotonous 
differential systems extensively studied by Hirsch et al. (see e.g. |12j). The behaviour of such 
competitive systems is known to be the most challenging, even when the equilibrium point is 
unique (which is not the case here). 

Both procedures are compared in the final section, using simulated and real data. Further 
numerical tests and applications are ongoing works in CA Cheuvreux. 

The paper is organized as follows: in Section [21 we make precise the modeling of splitting orders 
among several venues in the framework of Dark pools, first in static then in a dynamic way. This 
leads to an optimization problem under constraints. In Section [3j we study the execution function 
of one dark pool and introduce the recursive stochasic algorithm resulting from the optimization 
problem. In Section U] we analyze in depth this algorithm (a.s. convergence and weak rate) when 
the "innovations" (data related to the orders, the executed quantities and the market price) are 
assumed i.i.d. In Section [5] we extend the a.s. convergence result to a more realistic framework 
where these innovations are supposed to share some appropriate averaging properties (e.g. satisfied 
by a-mixing processes satisfying Ibragimov's condition). Section[6]is devoted to the second learning 
procedure, based this time on reinforcement principle, introduced in [3]. We make a connexion with 
the theory of (competitive) monotonous dynamical systems. Finally, in Sectional we present several 
simulations results on simulated and real data to evaluate the performances of both procedures with 
respect to an "oracle" strategy of an "insider" who could know a priori the executed quantities by 
every dark pool. 

Notations: • For every N > 1, set 1 N := {1,2,..., N}, V N := {r = (r;)i<;<„ G R+ | n = 

1}. Let l x := {ue R N | £ iG x «i = 0}. 

• Sij denotes the Kronecker symbol. 

• ( .| .) denotes the canonical inner product on M, d and | . | the derived Euclidean norm. 

• int(A) denotes the interior of a subset A of R d . 

• 5 a denotes the Dirac mass at a€ R . 
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2 A simple model for the execution of orders by dark pools 



2.1 Static modelling 

As mentioned in the introduction, we will focus in this paper on the splitting order problem in the 
case of (competing) dark pools. The execution policy of a dark pool differs from a primary market: 
thus a dark pool proposes bid/ask prices with no guarantee of executed quantity at the occasion 
of an over the counter transaction. Usually its bid price is lower than the bid price offered on the 
regular market (and the ask price is higher). Let us temporarily focus on a buying order sent to 
several dark pools. One can model the impact of the existence of N dark pools (N > 2) on a 
given transaction as follows: let V > be the random volume to be executed and let 9i £ (0, 1) be 
the discount factor proposed by the dark pool iG {1, . . . , N}. We will make the assumption that 
this discount factor is deterministic or at least known prior to the execution. Let ri denote the 
percentage of V sent to the dark pool i for execution and let D{ > be the quantity of securities 
that can be delivered (or made available) by the dark pool i at price OiS where S denotes the bid 
price on the primary market (this is clearly an approximation since on the primary market, the 
order will be decomposed into slices executed at higher and higher prices following the order book). 
The rest of the order has to be executed on the primary market, at price S. Then the cost C of 
the executed order is given by 



where pi = 1 — &i > 0, i = 1, . . . , N. At this stage, one may wish to minimize the mean execution 
cost C, given the price S. This amounts to solving the following (conditional) maximization problem 



However, none of the agents being insiders, they do not know the price S when the agent decides 
to buy the security and when the dark pools answer to their request. This means that one may 
assume that (V, D±, . . . , D n ) and S are independent so that the maximization problem finally reads 



where we assume that all the random variables min(V, D±), . . . , min(V, D N ) are integrable (other- 
wise the problem is meaningless). 

An alternative choice could be to include the price S of the security into the optimization which 
leads to the mean maximization problem 




i=i \ i=l 

/ N \ 





(2.1) 




(2.2) 




(2.3) 



4 



(with the appropriate integrability assumption to make the problem consistent). It is then conve- 
nient to include the price S into both random variables V and Di by considering V := V S and 
Di := DiS instead of V and Di which leads again to the maximization problem (|2.2|) . 

If one considers symmetrically a selling order to be executed, the dark pool is supposed to 
propose a higher ask price 9iS, 9i > 1, than the order book. The seller aims at maximizing the 
execution global (mean) price of the transaction. This yields to the same formal optimization 
problem, this time with pi = 9{ — 1, i = 1, . . . , N. 

All these considerations lead us to focus on the abstract optimal allocation problem (|2.2p which 
explains why the price variable S will no longer appear explicitly in what follows. 



2.2 The dynamical aspect 

In practice, there is no a priori assumption - or information available - on the joint distribution of 
(V, D±, . . . , D N ) under P. So the only reasonable way to provide a procedure to solve this allocation 
problem is to devise an on-line learning algorithm based on historical data, namely the results of 
former transactions with the dark pools on this security executed in the past. This underlines that 
our agent dealing with the dark pools is a financial institution like a bank, a broker or possibly a 
large investor which often - that means at least daily - faces some large scale execution problems 
on the same securities. 

This means that we will have to make some assumptions on the dynamics of these transactions 
i.e. on the data input sequence (V n , D™, . . . , -D™ )«>i supposed to be defined on the same probability 
space (Q,A,F). 

Our basic assumption on the sequence (Df, V n , i = 1, <, N) n >\ is of statistical - or ergodic - 
nature: we ask this sequence to be z^-averaging (a.s. and in L P (P)), at least on bounded continuous 

^ +1 ,fior(Mf 

assumption: 



functions, where v is a distribution on (M^ +1 , Bor(W+ +1 )). This leads to the following formal 



(i) the sequence (V n , D™, i = 1, . . . , N) n >i is averaging i.e. 



(ERG) U = { 



1 

k=l 



(ii) sup n E(V n ) 2 < +oo 



piV+l-, 



where =^ denotes the weak convergence of probability measures on IR + . For convenience, 
we will denote (V, Di, . . . ,D N ) the canonical random vector on IR^ +1 so that we can write v = 
£(y,D u ...,D N ). 

Assumption (ii) on the marginal distribution of the sequence (V n ) n >i is mainly technical. In 
fact standard arguments from weak convergence theory show that combining (£) and (ii) implies 

1 n 

-Vl/ fc ^Ey as n^oo 
n ^ 



(sup n K(V n ) 1+£ < +oo would be enough). An important subcase is the the (IID) setting 

(i) the sequence (V n ,Di, . . . ,Z?™) n >i is i.i.d. with distribution v = C(V,D\, . . . ,D N ), 

(ii) VeL 2 (F). 



(IID) 
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This more restrictive assumption is undoubtedly less realistic from a modeling point of view 
but it remains acceptable as a first approximation. It is the most common framework to apply 
the standard Stochastic Approximation machinery (a.s. convergence, asymptotically normal fluc- 
tuations, etc). So, its interest may be considered at least as pedagogical. The (ERG) setting is 
slightly more demanding in terms of assumptions and needs more specific methods of proof. It 
will be investigated as a second step, using some recent results established in [18] which are well 
suited to the specificities of our problem (in particular we will not need to assume the existence of 
a solution to the Poisson equation related to the procedure like in the reference book [3]). 

3 Optimal allocation: a stochastic Lagrangian algorithm 
3.1 The mean execution function of a dark pool 

In view of the modeling section, we need to briefly describe the precise behaviour of the mean 
execution function tp : [0, 1] — > P + of a single dark pool. 

Let (V, D) be an Revalued random vector defined on a probability space (O,, A, P) representing 
the global volume to be executed and the deliverable quantity (by the dark pool) respectively. 
Throughout this paper we will assume the following consistency assumption 

V > P-a.s. and ¥(D > 0) > 0. (3.1) 

The a.s. positivity of V means that we only consider true orders. The fact that D is not 
identically means that the dark pool does exist in practice. The "rebate" coefficient p is specific 
to the dark pool. 

To define in a consistent way the mean execution function of a dark pool we only need to 
assume that V G L X (P) (although more stringent integrability assumptions are made throughout 
the paper). 

Here the mean execution function ip : [0, 1] —> IR+ of the dark pool is defined by 

Vr€[0,l], (p(r) = pE(mm(rV, D)) (3.2) 

where p > 0. The function ip is finite, non-identically 0. It is clearly a concave non-decreasing 
bounded function. Furthermore, one easily checks that its right and left derivatives are given at 
every r 6 [0, 1] by 

ifi(r) = pE (l{r-v<D}^) and tp' r (r) = pE (l {rV<D} V) . (3.3) 

In particular, 

¥>'(()) = P E(V1 [D>0} ) > 

and if 

the (right continuous) distribution function of y is continuous on R+, (3-4) 

then 

(p is everywhere differentiate on the unit interval [0, 1] with ip' = ip\ on (0, 1]. 

Assumption (|3.4|) means that the distribution of y has no atom except possibly at 0. It can be 
interpreted as the fact that a dark pool has no "quantized" answer to an order. 

More general models of execution functions in which the rebate p and the deliverable quantity 
D may depend upon the quantity to be executed rV are briefly discussed further on. 
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3.2 Design of the stochastic Lagrangian algorithm 

Let V be the quantity to be executed by N dark pools. For every dark pool i £ X N the available 
quantity Di is defined on the same probability space (Q,A,F) as V. We assume that all couples 
(V,Di) satisfy the consistency assumption (|3.ip . 

To each dark pool i 6 I N is attached a (bounded concave) mean execution function ipi of 
type (j3.2|) . introduced in Section [27T] or (|8.ip . (|8.3[) studied in Section [8l 

Then for every r = (rj , . . . , r N ) € "P^ , 

AT 

SKn,...,^):^^). (3.5) 

i=l 

In order to design the algorithm we will need to extend the mean execution function <p (whatever 
its form is) as a concave function on the whole real line by setting 

(p(r) = - y^'(0) ifr<0 and tp(r) = <p(l) + logr if r > 1. (3.6) 

Based on the extension of the functions (pi defined by (|3.6|) . we can formally extend $ on the 
whole affine hyperplane spanned by V N i.e. 

U N := {reR N \^ = l}- 

i 

As announced, we aim at solving the following maximization problem 

max 3>(r) 

reV 

N 

but we will also have to deal for algorithmic purpose with the same maximization problem when r 
runs over % N . 

Before stating a rigorous result, let us a have a look at a Lagrangian approach that only takes 
into account the affine constraint that is max$(r) — A r{. Straightforward formal computations 

i 

suggest that 

r* G argmax-p^<I> iff <^(r*) is constant when i runs over X N 

or equivalently if 

1 r 

Vi£l N , tpiffi) = jr^2<p'j(rj). (3.7) 

3=1 

In fact this statement is not correct in full generality because the Lagrangian method does not 
provide a necessary and sufficient condition for a point to be a maximum of a (concave) function; 
thus, it does not take into account the case where $ reaches its maximum on the boundary &P N 
where the above condition on the derivatives may fail. So, an additional assumption is necessary 
to make it true as established in the proposition below. 
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Proposition 3.1 Assume that (V,Di) satisfies the consistency assumptions \3.1\) and (fff.^[ ) for 

every i £ I N . 

(a) Assume that the functions ip± defined by Ii3.2\) satisfy the following assumption 

Then argmax-p^ is a compact convex set and 

argmaxp^ = {re V N , \ tp'^n) = ^i( r i), i = 1, • ■ • ,N}. 
Furthermore argmax-^^$ = argmax-p^ . 

(b) If the functions ipi satisfy the slightly more stringent assumption, 



(C<) = min ip -(0) > max ^ ( — — - 



argmax-^<3? = argmax-p^ C int(V N ). 

Remarks. • If N = 2, one checks that Assumption (C) is also necessary to derive the conclusion 
of item (a). 

• As a by-product of the proof below we have the following more precise result on the optimal 
allocation r*: if r* € argmax-p^ and Ioi r *) '■= {} £ ^ N I r \ — 0}, then 

.max ^(0) < min <^(0). 

teZ (r*) t£Zo(r *) c 



Interpretation and comments: • In the case of a "regular" mean execution function, Assump- 
tion (C) is a kind of homogeneity assumption on the rebates made by the involved dark pools. If we 
assume that P(-Dj = 0) = for every i G I N (all dark pools buy or sell at least one security with 
the announced rebate), then (C) reads 

. > / Evl i^<pA 

EV ) 

since y^(0) = piKV. In particular, 

Assumption (C) is always satisfied when all the pi's are equal 
(all dark pools propose the same rebates). 

• Assumption (C) is in fact our main assumption in terms of modeling. It may look somewhat 
difficult to satisfy when the rebates are not equal. But the crucial fact in order to preserve the 
generality of what follows is that it contains no assumption about the dependence between the volume 
V and the "answers" Di from the dark pools. 
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Proof, (a) The function $ is continuous on a compact set hence argmax-p^ is not empty. Let 
r G argmax-p^<!> and Zo(r) := {i£ I N \ r{ = 0}. Clearly Zo(r) / I N so that card 2o(r) < N — 1. Let 
«6 l 1 such that Ui > 0, z G Xo(r). Then t i-)- <£(r + tu) defined on the right neighbourhood of 
reaches its maximum at so that its derivative at is non-positive. Specifying the vector u yields 

Viel (r), VjGX (r) c , ^(0) < ^(rj). 

Now if uG I -1- with U{ = 0, iG Zo( r )> then the t i— > $(r + tu) is defined on a neighbourhood of 
and reaches its maximum at t = so that its derivative is at 0; specifying the vector u yields 

Vi, je X (r) c , ^(r<) = ^-(rj)- 
Now, there exists at least one index i\ G Xo(r) c such that > |7q7^i > ]v^rT- Hence (^(r^) < 

^(jv^t) which implies in turn that for every i G 2b (r), ^ (0) < </4( r «) < V^dv^i)- Finally 
Assumption (C) implies that these inequalities hold as equalities so that 

VieJ N , ^(r-i) = y>i(ri). 

Conversely, let r G 7^ satisfying the above equalities. Then, for every r' G V N , the function 
1 1 — y <&(tr' + (1 — t)r) is concave on [0, 1] with a right derivative equal to at t = 0. So it is maximum 
at t = i.e. $(r) > $(r')- 

Now we pass to the the maximization over T~L N . Since it is an affine space and <3? is concave, it 
is clear, e.g. by considering $ as a function of {r\, . . . ,1"jv_i) 3 that 

argmax^* = {rG 7^, | tp'^n) = </?'i(n), i = l,...,N} 

(which is non-empty since it contains at least argmaxp^ ). Now let r G % N \V N . Assume there exists 

io G X N such tl 
Consequently 



1— V 1 

io G such that rj < 0. Then there always exists an index i\ G I N such that > ^J 8 " > j^zf- 



Vi (r i0 ) = (1 - r i0 )^ o (0) > ^ Q (0) > min^(0) > max^ f^y) > < (jyTTl) ^ 

which contradicts the equality of these two derivatives. Consequently all r^s are non-negative so 
that rG TV 

(6) If C< holds, the above proof shows that Xo(r) = so that argmax-p ^ C int(V N ). □ 



3.3 Design of the stochastic algorithm 

Now we are in position to devise the stochastic algorithm for the optimal allocation among the 
dark pools, taking advantage of the characterization of argmax-p <!>. In fact we will simply use 
the obvious remark that N numbers a\ v . . , a N are equal if and only if they are all equal to their 
arithmetic mean ° 1+ +ajv 



N 

We consider the mean execution function as defined by (|3.2[) . We assume from now on that the 
continuity assumption (|3.4p holds so that the representation (|3.3p of its derivative can be taken as 
its right or its left derivative on (0, 1] (and its right derivative only at 0). 
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Using this representation (|3.3p for all the derivatives ip\ yields that, if Assumption (C) is satisfied, 



then argmax-^^$ = argmax^$ and 



r*e argmax Pjv $ Vie {l,...,N}, E UHl{ r rv<Di} ~ JqYj , i' L {r'tV<D j } 

However, the set V N is not stable for the "naive" zero search algorithm naturally derived from the 
above characterization, we are led to devise the procedure on the hyperplane H N . 
Consequently, this leads to devise the following zero search procedure 

r n = r n - 1 + ln H(r n -\V n ,D?,...,D n N ), n > 1, r°eV N , (3.8) 

where, for every i£ I N , every r £ H N , every V > and every D\, . . . , D N > 0, 

1 N 

H i (r,V,D 1 , . . . ,D N ) = y^pdinV^D^nin&lo,!}} ~ j^^2Pj 1 {r J v<D j }n{r J e[o,i}} (3-9) 



3=1 

+R i (r,V,D 1 ,...,D N f) 

and the "innovation" (V n , D", . . . , -D™)n>i is a sequence of random vectors with non negative 

components such that, for every n > 1, (V n ,Df ,i = 1,<,N) = (V,Di,i = 1,<,N) and the 
remainder terms Ri have a mean-reverting effect to pull back the algorithm into V N . They are 
designed from the extension (|3.6|) of the derivative functions <*p\ outside the unit interval [0, 1]; to 
be precise, for every i£l N , 

Ri(r, V,D 1 ,...,D N ) = pi ^(1 - r^l^.jojnjr^o} + — IjvxDJn-^i}^ 

1 N ( 1 A 

3.4 Interpretation and implementability of the procedure 

[> Implementability. The vector {r2)i<i<N in (|3.8p represents the dispatching of the orders 
among the N dark pools to be sent at time n + 1 by the investor. It is computed at time n. On the 
other hand V n represents the volume to be executed (or its monetary value if one keeps in mind 
that we "plugged" the price into the volume) and the D™ the "answer" of dark pool i, still at time 
n. 

The point is that the investor does have no access to the quantities Df. However, he/she knows 
what he/she receives from dark pool i, i.e. min(Z)™, r™ _1 l/ n ). As a consequence, the investor has 
access to the event 

{minp^r™- 1 ^) = r^V" 1 } = {r^V™ < Df} 

which in turn makes possible the updating of the procedure although he/she has no access to the 
true value of Df . 
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So, except for edge effects outside the simplex V N , the procedure as set can be implemented on 
real data. 

> Interpretation. As long as r is a true allocation vector, i.e. lies in the simplex V N , the 
interpretation of the procedure is the following: assume first that all the factors pi are equal (to 1). 
Then the dark pools which fully executed the sent orders (r»V < Dj) are rewarded proportionally to 
the numbers of dark pools which did not fully executed the request they received. Symmetrically, the 
dark pools which could not execute the whole request are penalized proportionally to the number 
of dark pools which satisfied the request. 

Thus, if only one dark pool, say dark pool 1, fully executes the request at time n, its pourcentage 
will be increased for the request at time n + 1 by 7„(1 — jf)V n i-e. it will asked to execute 
r" = r™ _1 + 7„(1 — jq)V n % of the total order V n+l . The other N — 1 dark pools will be penalized 
symmetrically: the pourcentage rf of the total request V n+1 each of them will receive at time n + 1 
will be reduced by "f n jfV n - 

If k dark pools totally execute their request at time n and the N — k other fail, the pourcentages 
of V n+1 that the "successful" dark pools will receive for execution at time n + 1 will be increased 
by 7„(1 — jf)V n , each of the N — k "failing dark pools" being reduced by r ) n j^V n . 

If no dark pool was able to satisfy their received request at time n, none will be penalized and 
if all dark pools fully execute the received orders, none will be rewarded. 

In short, the dark pools are rewarded or penalized by comparing their mutual performances. 
When the "attractivity" coefficents pi are not equal, the reasoning is the same but weighted by 
these attractivities. 

[> Practical implementation. One may force the above procedure to stay in the simplex V N 
by projecting, once updated, the procedure on V N each time it exits the simplex. This amounts to 
replace the possibly negative rj by 0, the rj > 1 by 1 and to renormalize the vector r by dividing 
it by the sum of its terms. 

Furhermore, to avoid that the algorithm leaves too often the simplex, one may simply normalize 
the step 7„ by considering the predictable step 

~ _ n - 1 _ Jn_ 

ln ~ ln X V 1 + • • • + V n ~ l ~ EV 

4 The (IID) setting: a.s convergence and CUT 

Theorem 4.1 Assume that (V,D) satisfy \3.1\) . that V £ L 2 (P) and that Assumption (C) holds. 
Assume furthermore that the distribution of y satisfies the continuity Assumption fl3.^[ ). Let 7 := 
(7n)n>i be a step sequence satisfying the usual decreasing step assumption 

^7n = +oo and ^ In < +°°- 

n>l n>l 

Let (V n , D™, . . . , D™ ) n >i be an i.d.d. sequence defined on a probability space (Q, A, P). Then, there 
exists an argmax-p^$ -valued random variable r* such that 

r n — > r* a.s. 

Lf the functions ipi satisfy (C<) then argmax-p $ C int(P N ). 
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Proof of the theorem. In this setting, the algorithm is (non homogenous) Markov discrete time 
process with respect to the natural filtration T n := <r(r°, {V k ,D\, . . . ,D^), 1 < k < n) with the 
following canonical representation 



r n+1 = r n + ln+1 H{r n ,V n+ \D n l + \...,D n N +l ), r° e V 



N 



= r n + 7n+i/i(r n ) + 7 n+ iAM n+ i 

where, for every r£ H N , 

h{r) := E H(r, V, D x , . . . , D N ) = L^n) - ^ 
is the so-called mean function of the algorithm, and 




Ki<N 



AM n+1 = H(r n ,V n+1 ,D^ + \...,D n N +1 )-E(H(r n ,V n+1 ,D^ + \...,D n N +1 )\T n ) 
= H{r n ,V n+1 ,D^ + \...,D n N +1 )-h(r n ) 

since (V n+1 ,D^ +1 , D n N +1 ) is independent of T n . 

One derives from Proposition I3,ll fa) that the mean function h of the algorithm satisfies {h = 
0} = argmax-p^ and that, for every r £ H N \{h = 0} and every r* £ {h = 0}, 

N 

(h(r) \r-r*) = (h(r) - h(r*) \ r - r*) = V (^(r<) - ^(r*))(r, - r*) < (4.1) 

simply because each function 99- is non-increasing which implies that each term of the sum is non- 
positive. The sum is not zero otherwise f'(ri) = ^(r?) as soon as ^ r* which would imply 
h(r) = 0. 

The random vector V being square integrable, it is clear that Hi(r, V,D±,..., D N ) satisfies the 
linear growth assumption 

Viel N , VreH N , \\H l (r,V,D 1 ,...,D N )\\ 2 <2(m a x Pj )\\V\\ 2 (N+\r\) 

3 

At this stage one may conclude using a simple variant of the standard Robbins-Monro Theorem 
(like that established in |20|): there exists a random variable r* taking values in {h = 0} such that 

c 



4.1 Rate of convergence 

Our aim in this section is to show that the assumptions of the regular Central Limit Theorem 
(CLT) for stochastic approximation procedures are fulfilled. For a precise statement, we refer 
(among others) to [3] (Theorem 13 p. 332). For the sake of simplicity, we will assume that the mean 
function h has a single zero denoted r*. The following lemma provides a simple criterion to ensure 
this uniqueness. 
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Lemma 4.1 Assume that all the functions ipi, iG T N , are decreasing (strictly). Then 

{h = 0} = argmaxp $ = r* G int(V N ). 

Proof. In particular (C<) is satisfied so that argmax-p^$ C r* G int(V N ). If r, r'G {fo = 0}, r 7^ r', 
it follows from (|4.ip that (p[{ri) = </^(rQ for some index i such that Ti^r\. □ 

The second ingredient needed to establish a CLT will be the Hessian of function To ensure 
its existence we will make one further assumption on a generic random couple (V, D), keeping in 
mind that P(D > 0) > 0, but that P(D = 0) may possibly be positive too. Namely, assume that the 
distribution function of (V, D) given {D > 0} is absolutely continuous with a probability density 
/ defined on (0, +00) 2 . Furthermore we make the following assmptions on /: 

(i) for every v > 0, u h-> f(v,u) is continuous and positive on (0,oo), 

o , (4-2) 
(it) Ve€(0,l), sup f(y,u)V 2 £L l {¥). 

eV<u<V/e 

Note that (ii) is clearly always satisfied when V G L 2 (P) and f is bounded. The conditional 
distribution function of D given {D > 0} and V is given by 

F D (u I V = v, D > 0) := ¥(D <u\V = v,D > 0) = / /(«, u')**', u > 0, u > 0, 

Jo 



Lemma 4.2 (a) Assume (V, D) satisfies the above assumption j^.i?| ). T/ien £/te mean execution 
function <p(u) := pE(min(iiV, D)) is concave, twice differentiable on P+ and for every u > 0, 

<p"(u) = - P E(V 2 l {D>0} f(V,uV)) < 0. 

(b) If (V,Di) satisfies the above assumption \4-*fy f or every iG I N , then the function h defined on 

by h(u\, . . . , u N ) = ( <p'i(ui) — if J2i<j<N vU^i) ) ^s differentiable on (0, oo)^ and admits 

a continuous extension on given by 

with ai(u) = —(p'-(u) > 0. 

l<i,j<N 

(c) Lei A := [— a,j + iVoi5y]i<ij<jv'> a i 5 ■ ■ ■ ,a- N > and Ze£ a = minjaj. Its kernel Ker(^4) is one 
dimensional, A(M> ) = l -1- andyl^ is bijective. Every non-zero eigenvalue A with eigenspace E\) 
satisfies 

R(X)>Nxa and £ A C l 1 . 

Proof, (a) is a straightforward application of the Lebesgue differentiation Theorem for expectation. 

(b) is a consequence of (a). 

(c) The transpose A* of ^4 has a strict dominating diagonal structure i.e. A 1 - > 0, A|- < 0, i ^ j 
and ^ ■ A^- = for every i. Consequently, it follows from Gershgorin's Lemma (see [9]) that is an 
eigenvalue of order 1 of A* (with 1 as an eigenvector and that all other eigenvalues have (strictly) 



Dh(u) = -- 



a,j(uj) + Nai(ui) 5 i: 
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positive real parts). Consequently Ker(A) is one dimensional. The fact that A(M. ) C 1^ is obvious 
so that this inclusion holds as an equality by the dimension formula. Hence all the eigenvectors 
not in Ker(j4) are in 1^. Set 5j — dj — a > 0, i = 1, . . . , N. Then A t has a dominating diagonal 
structure so that all its eigenvalues have non-negative real parts. Now if A is an eigenvalue of A, it 
is obvious that A — Na is an eigenvalue of A. Consequently 3ft(A) > Na. □ 



Theorem 4.2 Assume that the assumptions of Theorem \4-l\ holds and that argmax<3? is reduced 
to a single point r* £ V N so that r n — > r* W-a.s. as n — > oo. Furthermore, suppose that Assump- 
tion gjfy holds for every (V, A), and that V G L 2+5 (F), 5 > 0. Set 

c 1 
7 n = — , n > 1 lOTWl c > 



' " 25fte(A min ) 

where A m i n denotes the eigenvalue of A 00 := — Dfo(r*).^_L uni/i i/ie lowest real part. Then 

' "A^O;!] 00 ) 



V7n 

where the asymptotic covariance matrix S°° is given by 

r-CO 

e u{A°°-g) c oc e u(A°°-gy du 







where 



C°° = E (H(r*, V, A, ■ ■ ■ , D N )H(r*, V,D 1 ,..., Dj*)^ 

and (A°° - |^)* stands for the transpose operator of A°° — ^ e £(1^). 

Remark. The above claim is consistent since u i— >■ iJ(r, v,5±, . . . ,5 N Y u preserves 1^. 

Proof. First note that, since r* 6 int(V N ), the above Lemma l4.2l ffr) shows that (still making the 
confusion between the linear operator Dh(r*) and its matrix representation in the canonical basis) 

Dh(r*) = -1 [-^(r*) + iVa;(r*) ... Y with a,(r) = ft E(F 2 l {A>0} /(F, rV)) > 

Then, Lemma 14.21 (c) implies that — Dh(r*)^± has eigenvalues with positive real parts, all lower 
bounded by minj cij(r*) > 0. 

At this stage, one can apply the CLT for stochastic algorithms defined on l -1- (see e.g. [3], 
P-341). n 

5 The (ERG) setting: convergence 

For the sake of simplicity, although it is not really necessary, we will assume throughout this section 
that 

argmax-p^ = {r*} C mt(7 7 JV ) 
possibly because all the execution functions <pi are decreasing so that, following the former Lemma l4.ll 
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So we assume that the sequence (V n ,D^,i = 1, . . . , N) n >i satisfies (ERG) U with a limiting 
distribution v such that, for every i £ Z N , its marginal v% = £(V, Di) satisfies the consistency 
assumption (|3.ip and the continuity assumption (|3.4|) . We will also need to make a specific as- 
sumption: there exists £q > such that 



P(V > eo) > 



(ii) suppf^Cf^, i = 1, . . . , N | {V > £q} ) ) is a neighbourhood of V N in 



(5.1) 



This assumption means that all allocations across the pools lying in the neihbourhood of V N can 
be executed. 

On the other hand, it follows from (ERG) U and some standard weak convergence arguments 
that 

1 n 2 

VieX^, VueM+, -^y fc l { ^ fe < Dn -E(yi {ny < A} ) a -^> L2 as n oo, 

fe=i 

since the (non-negative) functions f u (v,5) := ulj u „<ai, u > 0, are P/^^-o.s. continuous and 0(i>) 
as v — > +oo by (j3.4h . Moreover this a.s. convergence holds uniformly on compact sets with respect 
to u since u \-¥ E,Vl{ u y < D i } is continuous, still owing to (|3.4p . Our specific assumption is to require 
a rate in the above a.s. and L 2 (P)-convergence. Namely, we assume that there exists an exponent 
cti £ (0, 1] such that 

1 n 

VneM+, -Y J V k l {uVk<D , } -nV\uV<D 1 }) = 0{n~^) a.s. and in L 2 (P). (5.2) 

k=l 

This assumption e.g. from the more general assumption that, for every i E In, the marginal 
Ui = C(V,Di) satisfies (|3.4p and 

(y n ,D™) is i/j-averaging at rate 

on a subspace V Qii 2 containing all the functions f u . 

Note that, when the sequence (V n ,Df,i = l,...,iV) n >i is i.i.d. with distribution v then 
elementary martingale arguments show that the whole sequence is ^-averaging at rate \ — n for 
every rj € (0, 1/2) on Vi 2 = ^ 2 ( y ) ( and a11 /« £ L 2 (v), u> 0, since V 6 L 2 (P)). So, the theorem 
below almost embodies the a.s. convergence theorem established in the (I ID) setting (except for 
the integrability assumption on V). 

Now we are in position to state the main convergence result of this section. We rely on the 
extension of Robbins-Siegmund Lemma proposed in |18j . For the reader's convenience it is recalled 
in the Appendix. 

Theorem 5.1 Let (V n , D™, . . . , D r ^) n >Q be a sequence of input satisfying {ERG) V and such that, 
for every i£ X N , the marginal distribution = C(V, Di) satisfies the consistency assumptions L3. 1\) 
and \3.4\ )- Suppose furthermore that, the sequence (V n , -Df)n>i satisfies the rate assumption i5.2\) . 
If the step sequence ( r y n ) n >i satisfies 

= +oo, 7 n = o{n-~ l ) and ^re 1- - max(7 2 , \y n - 7 n+ i|) < oo 

n>l n>l 

where a := muijgx^ on G (0,1], then the algorithm defined by §3. 9\) a.s. converges towards r* = 
argmaxp . 
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technical comment. The above condition on the step sequence (7 n )n>i is satisfied as soon as 
In = with /3g (I - a, 1]. 

Proof. Step 1. First, we aim at applying the extended Robbins-Siegmund Lemma established 
in [18] (see Appendix, Theorem A.l for its statement) for stochastic algorithms with ^-averaging 
inputs dynamics in presence of a Lyapunov function. We will consider the case p = 2 and /3G (0, a]. 
We set G = -H and AM" = and we consider the input Y n = (V n+1 , D™ +1 , D n N +1 ), n > 0. 
Let L(r) = ||r — r*| 2 be our candidate as a Lyapunov function. 

First note that it follows from (13. 9ft that the function H satisfies the growth assumption flA.4|) 
since 

Vr G n n , Vy G M JV+1 , |ff(r,y)| < C^(y)(l + |r|) 

where Ch > and g(v, 6%, . . . , 6 N ) = v. 

In view of the ergodic assumption (|5.2p and the fact that r* lies in V N , it is clear from its 
definition that H(r*,.)£ Vp t 2 for every /3G (0,a]. 

At this stage it remains to check the "weak local Lyapunov" assumption ()A.5p for G = —H. 
This fact is obvious since, for every r G 7^ and every input y = (y, Si, . . . , 5n) G (0, +oo) x M. N , 

N 

(H(r,y) - H{r\y)\r - r*) = - H^r*, v, «5 4 ))(r 4 - r*) < 

i=i 

where 

Hi(u,V,5i) = PiV fl^^l^l] (") + (!- u)l Si>O ,u<0 + ^-v<Si,u>l) ,i£l N (5.3) 

is clearly non-increasing with respect to u. 

At this stage, using that sup n>1 E(V n ) 2 < +00, we can apply our extended Robbins-Siegmund 
lemma also that 



> Loo < +00 a.s. and ^ j n (r n - r* | G(r n , Y n ) - G(r*, Y n )) < +00 a.s. (5.4) 



n>l 



Step 2. At this stage it suffices to show that r* is a.s. a limiting point of (r n ) n >o since \r n — r*\ 
converges to < +00 a.s. 

Let 77 denote a positive real number such that, for every i£ I N , [r* — r],r* + rj\ C (0, 1). One 
derives from (|5.3p and the monotonicity of Hi(u,v,5) in ti£ K that for every i G X N and every 
r€H N , 

(Hi(ri,v,5i) - Hi{r*,v,Si))(r* - n) > Pivr]l {ri>r * +v} l {s/veJv} 
where J v = (r*,r* + rj). As a consequence 

(G{r,y) - G(r*,y) \r-r*)> E £^l{e> eo }l W 60, J^l ri >rj+»j- 

where p = min^ pi and the open set 0„ is defined by 

O v = jy = . ,<5jv)G (eo,+oo) x s.t. ^ £ J„, i £ X^j . 
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Now, one derives from (|5.4p that 

^2lnl0 v (Y n ) 1 r^>r*+r ) < +OO a.S. 

n iei N 

Now Assumption f)5. 1 1) implies that v{0^) > 0. Furthermore v(dO v ) = owing to the continuity 
assumption so that (ERG) U implies 

1 n 

-Ylo v (Y k ) — ► V (0„) > a.s. 
fe=l 

An Abel transform and the facts that the sequence 7 n is non-increasing and Yln>i In = +°o 
classically implies that 

^7«(lo, ; ( yfc ) - v(O n )) a.s. converge 

n>l 

so that 

^ 7n l 0l) (y n ) = +oo a.s. 

n 

In turn, this implies that 

liminf l{ r n >r * +7? j = a.s. 

This holds of course for a sequence of real numbers rf decreasing to 0. 

Let 7£oo be the set of limiting values of the sequence (r n ) n >o. It is is a.s. non-empty since the 
sequence (r n ) ra >o is bounded. Then T^oo is a.s. compact and it follows from what precedes that 
T^oo H rii<i<7v ( — 00 ) r i + V i ¥^ ® (and is compact). Hence a decreasing intersection of non-empty 
compact sets being a (non-empty) compact set T^oo n rii<i<Ar( — °°> r i\ / On the other hand 
T^oo C 7i N since the algorithm is "H^-valued. But ni<i<Ar( — °°; r i \ ^^jv = i r *}- Consequently r* 
is a limiting point of the algorithm which implies that it is its true a.s. limit. □ 

Application to o-mixing stationary data. If (V n ,Df,i = l,...,N) n >i is a stationary a- 
mixing sequence which mixing coefficients (a n ) n >i satisfy Ibragimov's condition for some 5 > 0: 

n>l 

(which is satisfied in case of geometric a-mixing) then the sequence (V n , Df,i = 1, . . . , N) n >i is v- 
averaging where v is the stationary marginal distribution of the sequence (supposed to satisfy ()3.ip 
and (|3.4p ) at rate (3 for every /3e (0, 1/2). To be precise, L 2 (v) C V + 2 an d 

l 2+5 (^)c p| v A2 . 

0</3<± 

In particular, all the functions f u (v, 5) := vls uv <gy, u > 0, lie in every Vpp, < /? < |, so that the 
rate condition (|5.2p is satisfied. 
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As concerns the stationary assumption on the input data sequence, it can be considered as 
realistic if one think of execution objectives given on a daily basis. 



Example: An exponential discrete time Ornstein-Ulhenbeck model for (V n , D™, . . . , D r 



yn = v O e XS^ D n = ^Xf ^ i = \ f , , , } N t n > 1, 

where v°, . . . , d° N are positive real numbers and the sequence (A n ) n >i satisfies the linear 
auto-regressive dynamics 

X n+1 = m + AX n + B E n+ \ n > 1, 

withm g R N+1 , Ae M(N+1,N + 1,M), \\A\\ < 1, B£ M(N+1,M) with mnk(B) = N+1 (< M) 
and (H n ) n >i is an i.i.d. sequence of AA(0; /^^-distributed random variables. We assume that the 
sequence is stationary i.e. that the distribution of A 1 is the (Gaussian) invariant distribution 
(with covariance matrix C solution to the Lyapunov equation C — AC A 1 = BB l where * stands for 
transpose). Then (see [6], p. 99), the sequence (A n ) n >o is geometrically a-mixing and subsequently 
so is (V n ,D™, . . . , -D") n >i (with respect to its natural filtration). Furthermore, it is clear that its 
distribution v satisfies the dispersion assumption (|5.ip the process [X\ y X\ — Aq,z = 1, . . . , N) is 
a non-degenerate Gaussian distribution over since B has full rank N + 1. 

6 An alternative procedure based on a reinforcement principle. 

Recently, inspired by the discussion developed by Almgren and Harts in [2] about liquidity esti- 
mation, Berenstein and Lehalle devised a "smart routing" recursive procedure of requests to be 
executed by a pool of N dark pools (see [4]). This procedure is not based on the opimization of a 
potential function but on a intuitive reinforcement mechanism. Let If be the profit induced by the 
execution of the order sent to dark pool i at time n. The proportion rf of the global order V n+l 
to be sent to dark pool i for execution at time n + 1 is defined as proportional to this profit i.e. by 

I n 

VieZ 



N ' i r> m' 

^3 3 

The updating of the random vector I n is as follows 

Vn>0, Viel N , ir +1 =Jf + p i mm{r?V n+1 ,D?+ 1 ), lf = 0. 

The first equation models the idea of "reinforcement" since the proportion of orders sent for ex- 
ecution to dark pool i is proportional to the historical performances of this dark pool since the 
beginning of the procedure. 

The second equation describes in a standard way - like in the optimization algorithm - the way 
dark pools execute orders. 

Elementary computations show that the algorithm can be written directly in a recursive way 
in terms of a new vector valued variable 

jn 

X n = —, n>l, 
n 



since 



18 



This is a standard form a stochastic algorithm (with step j n 
Furthermore, note that setting p := rninj Oj, 



^/r>pmin(ly",minA n ) 

ltX JV 



since rf > for at least one dark pool iE X N . Consequently, as soon as the sequence (V n , D™, . . . , D" ) 
is stationary and ergodic 

1 JL /I . A 

— V, min IX- a.s 



limmf X? > ohm - Vmin ( —V , mm L> fc | = pEmin ( — V, 

iex k=i v jy / \ 



So if we make the natural assumption that 

Emin f — V. min A I > 

then, a.s., the algorithm X n cannot converge to 0. 

If we make the additional assumption that the sequence (V n , -D™, . . . , D™ ) is i.i.d. then the 
algorithm is a discrete time (non homogenous) J^-Markov process with respect to the filtration 
T n = o{V k ,D\, . . . , D* k = 1, . . . , n), n > 0, so that it admits the canonical representation 

x n+l =x n_ 7n+i(x n_^ (r n ))+7n+iAM n+l W > 0, 

where 7 n = - and 

AAf" = p J min(rrV ra ,A n )-^K"" 1 ), *G X w , n > 1, 

is an J-" n -martingale increment. Furthermore it is L 2 -bounded as soon as L 2 . 

In fact the specific difficulties induced by this algorithm are more in relation with its mean 
function 

h : x i — > [Xi-ipil X * (6.1) 

than with the martingale "disturbance term" r y n+ \A.M n+1 . Our first task will be to prove under 
natural assumptions the existence of a non degenerate equilibrium point. Then we will show why 
this induces the existence of many parasitic equilibrium points. 

6.1 Existence of an equilibrium 

In this section, we will need to introduce a new function associated to a generic order V and a 
generic dark pool with characteristics (p,D). 

4(u):=^, u>0, il>(0)=v'(0) = P EVl {D>0} . (6.2) 
u 

If Assumption (13. ip holds then ^>(0) < +oo and ip is continuous at 0. It follows from the concavity 
of <p and f(0) = that ^ is non-increasing. It is continuous as soon as <p is e.g. if Assumption (|3.4p 
holds true. 
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Proposition 6.1 Let N > 1. Assume that Assumption \3. 1\) holds for every couple (V,Di), iG T N . 

(a) There exists a x* G smc/i i/tcrf 

J>*>0 and Vil ^*' A =xh ieV (6-3) 
iei N \^j^N X iJ 

(b) Let ipi be the functions associated to dark pool «G I N by $6.2\) . Assume that for every i£ T N , 
ipi is (continuous and) decreasing on [0,oo) and that 

VT 1 (min^(0))<l. (6.4) 

Then there exists x* G int(V N ) satisfying 116. ,S\) . 
Proof, (a) We define for every r = (ri, . . . , r N ) G 7^ 



This function maps the compact convex set V N into itself. Furthermore it is continuous since, on 
the one hand, for every i£ I N , (fi is continuous owing to the fact that (V, Di) satisfies (|3.ip and, 
on the other hand, 



3&N 

Indeed, for every i£l N , 



J2 ^i( r i) ^ ™m <pj (j^j > 0. 



Vi(^) ^ E min(y,A)>0 



since V > P-a.s. and P(-Dj = 0) < 1. Then it follows from the Brouwer Theorem that the 
function has a fixed point r* . Set for every i£ I N , 

It follows immediately from this definition that 

VieT N , x* = Lpi{r*) 

which in turn implies that V^„ cT Viirf) = YlncT x *n so that r* = ^ Xi — «, iG T N . 

(6) For every i£ T N we consider the inverse of ipi defined on the interval (0, </^(0)]. This function is 
decreasing continuous and lim^o V^ -1 ^) = +oo. Then, let be the continuous function defined 
by 

V0G(O,mW(O)], e(9)=Y,^ 1 (0). 

' fcX JV 
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We know that limg^o ©(#) = +00 and we derive from Assumption (|6.4p that 6(mnij e z ¥>i(0)) < 1. 
So, owing to the (strict) monotonicity of 6, there exists 9* £ (0, min g x 99^ (0)) such that @(6*) = 1. 



Then r* := (rjf, . . . , r* ) € int{V N ) since £\ r| = 1 by definition of (9*. If r* Q = 0, then 9* = ip io (0) = 
minjgj^ t/^(0) which is impossible. □ 

Corollary 6.1 Assume that all the functions ipi are continuous and decreasing. If furthermore, 
the rebate coefficients pi are equal (to 1) and ifW(Di = 0) = for every i€. I N then there exists an 
equilibrium point lying in int(V N ). 

Proof. Under the above assumptions </^(0) = EF > 0. Consequently 



Comments. Unfortunately there is no hope to prove that all the equilibrium points lie in the 
interior of V N since one may always adopt an execution strategy which boycotts a given dark 
pool or, more generally, iVo dark pools. So it seems hopeless to get uniqueness of the equilibrium 
point. To be more precise, under the assumptions of claim (6) of the above Proposition 16.11 there 
exists at least one strategy involving a subset of N — Nq dark pools Nq = 0, . . . , N — 1 (one dark 
pool is needed at least). Elementary combinatorial arguments show that there are at least 2 N - 1 
equilibrium points. 

So, from a theoretical point of view, we are facing a situation where there may be many parasitic 
equilibrium points, some of them being clearly parasitic. However it is quite difficult to decide a 
priori, even if we make the unrealistic assumption that we know all the involved distributions, 
which equilibrium points are parasitic. 

This is a typical situation encountered when dealing with procedures devised from a reinforce- 
ment principle. 

However, one may reasonably hope that some of them are so-called "traps" , that means equi- 
librium points which are repulsive at least in one noisy direction so that the algorithm escapes from 
it. Another feature described below suggests that a theoretical study of the convergence behaviour 
of this procedure would need a specific extra work. 

The next natural question is to wonder whether an equilibrium x* of the algorithm - namely 
a zero of h - is (at least) a target for the algorithm i.e. is attractive for the companion ODE, 
x = —h(x). 

Proposition 6.2 An equilibrium x* satisfying $6.3\) is locally uniformly attractive as soon as 



Set 



C 1 (min^(0))=C 1 (^(0))=0<l. □ 



'JV 
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Remark. In fact the following inequalities are satisfied by any equilibrium x*: 

This follows from the convexity of the function £ i— >■ £ — PilStj which is zero at x* with positive 
derivative. As a consequence, the right hand side in the above sufficient condition is always positive 
which makes this criterion more realistic. 

Proof. Elementary computations show that the differential Dh(x) of h at x G is given by 

mez,, ^)-%(l-i»5(£)) + ^(f). 

As a consequence all the diagonal terms of Dh(x*) are positive. The above condition for all the 
eigenvalues of Dh(x) to have positive real parts follows from a standard application of Gershgorin's 
Lemma to the transpose of Dh{x). □ 

6.2 A competitive system 

But once again, even if we could show that all equilibrium points are noisy traps, the convergence 
would not follow for free since this algorithm is associated to a so-called competitive system. A 
competitive differential system x = h(x) is a system in which the field h : W N — > WL N is differentiable 
and satisfies 

VxgR n , Vi, j el N , i^j, ^%)>0. 

As concerns Almgren and Harts's algorithm, the mean function h is given by (|6.ip . and under 
the standard differentiability assumption on the functions ipi's, 

VxeR N , p-(x) = <f/ i ( ^ )- f ^ ^>0. 

These systems are known to have possibly a non converging behaviour even in presence of a 
single (attracting) equilibrium. This is to be compared to their cooperative counterparts (with 
negative non-diagonal partial derivatives) whose flow converge uniformly on compact sets toward 
the single equilibrium in that case. This property can be transferred to the stochastic procedure by 
the mean if the so-called ODE method which shows that the algorithm almost behaves like some 
trajectories of the Ordinary differential Equation associated to its mean field h (see e.g. \13\ EJ [7] 
for an introduction). Cooperativeness and competitiveness are in fact some criterions which ensure 
some generalized monotonicity properties on the flow of the ODE viewed as a function of its space 
variable. For some background on cooperative and competitive systems we refer to [ 1 1 1, 1 1 2 j and the 
references therein. 

7 Numerical tests 

The aim of this section is to compare the behaviour of both algorithms on different data sets : 
simulated i.i.d. data, simulated a-mixing data and (pseudo-)real data. 
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Two natural situations of interest can be considered a priori: abundance and shortage. By 
"abundance" we mean (in average, the requested volume is lower than the available 

one). The "shortage" setting is the reverse situation where EV > "Yld^EDi. 

In fact, in the "abundance" setting, both our procedures (optimization and reinforcement) tend 
to remain "frozen" at their starting allocation value (usually uniform allocation) and they do not 
provide a significant improvement with respect to more naive approaches. By contrast the shortage 
setting is by far more commonly encountered on true markets and turns out to be much more 
challenging for our allocation procedures, so from now on we will focus on this situation. 

Our first task is to define a reference strategy. To this end, we introduce an "oracle strategy" 
devised by an insider who knows all the values V n and D™ before making his/her optimal execution 
requests to the dark pools. It can be described as follows: assume for simplicity that the rebates 
are ordered i.e. p\ > p2 > ■ ■ ■ > Pn- Then, it is clear that the "oracle" startegy yields the following 
cost reduction (CR) of the execution at time n > 1, 



CR' 



oracle 




io-1 \ io— 1 io 

i=l / i=l i=l 

N 

if y, d ? < yn - 



i=l i=l 
Now, we introduce indexes to measure the performances of our recursive allocation procedures. 

• Relative cost reduction (w.r.t. the regular market): they are defined as the ratios 
between the cost reduction of the execution using dark pools and the cost resulting from an 
execution on the regular market for the three algorithms, i.e., for every n > 1, 

Q j^oracle 



o Oracle: 



CR algo J2ZlPi min (. r i Vn > D i 



o Recursive "on-line" algorithms: 
(with algo = opti, reinf). 



yn yr, 



Performances (w.r.t. the oracle): the ratios between the relative cost reductions of our 
allocation algorithms and that of the oracle, i.e. for every n > 1 

CR opti CR reinf 

and 



(J R oracle ' (J R oracle 

which seems a more realistic measure of the performance of our allocation procedures since the 
oracle strategy cannot be beaten. 

Since these relative cost reductions are strongly fluctuating (with variables V n and Df in fact), 
we will plot the moving average of these ratios (on the running period of interest) and express them 
in pourcentage. 

Moreover, when we simulate the data, we have chosen 10 4 simulations because it corresponds 
approximatively to the number of pseudo-real data observed within a day. 

The choice of the gain parameter is the following (in the different settings considered below) 



Q 

7n = -, n>l 
n 



where c equals to some units. 
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7.1 The (I ID) setting 



We consider here simulated data in the i.i.d. setting, where the quantity V and Dj, i £ X/v, are 
log-normal variables and N = 3. The variables V and Di, i £ In, satisfy the assumptions of the 
CLT and we have the rate of convergence at least of the optimization algorithm. 
The shortage setting is specified as follows: 



N 



i=l 



with 

/o.or 

EDi = i,l<i<N, Var(F) = l,Var(A) = M < * < N and p = 0.03 

\0.05, 

The running means of the performances are computed from the very beginning for the first 100 
data, and by a moving average on a window of 100 data. 

The initial value for both algorithms is set at r° = 1 < i < N . 



Relatives Cost Reductions 




Performances 



Oracle 

Optimization/Oracle 

- - - Reinlorcement/Oracle 



Number of simulations x 10 



10 20 30 40 50 60 70 
Number of simulations x 10 2 



Figure 1: Shortage setting Case N = 3, ray = | Yli=i m D i ^ m Di = h o~y = 1, au i = 1, 1 < i < N. 

As expected, the optimization procedure outperforms the reinforcement one and both proce- 
dures quickly converge (see Figure [1]) with respect to the data set size. Note that the allocation 
coefficients (not reproduced here) generated by the two algorithms are significantly different. A 
more interesting feature is that the performances of the optimization procedure almost replicate 
those of the "oracle" . Further simulations suggest that the optimization algorithm also seems more 
robust when the variances of the random variables fluctuate. 



24 



7.2 The (ERG) setting 



We consider here simulated data in the ergodic setting, where the quantity V and Dj, i£ In, are 
exponentials of an Ornstein-Uhlenbeck process, i.e. 



X n+l = m + AX n + BE n+ \ 



where \\A\\ < 1, BB* G GL(d,R) and 



to 




3JV+1 



\"iV+l/ 



~N(0,I N+1 ) i.i.d., e 



A' 



We are still interested in the shortage situation. The initial value of the algorithms is r. 
1 < i < N and we set 



— J_ 

N ' 



/0.01N 




(A 




0.03 , 


m = 




I— 


\0.05/ 
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The running means of the performances are computed from the very beginning for the first 100 
data, and by a moving average on a window of 100 data. 



Relative Cost Reductions 
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— - Reinforcement 
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Optimization/Oracle 

Reinforcement/Oracle 



20 40 60 80 

Number of simulations x10 2 



Figure 2: Shortage setting: Case N = 3, my > YliLi m D t , o~v = 1-21, (Xd = (8.21, 3.05, 1.07)'. 

We observe in this ergodic setting a very similar behaviour to the i.i.d. one, with maybe a more 
significant advantage for the optimization approach (see Figure [2] right): the difference between the 
performances of both algorithms reaches 11% in favour of the optimization algorithm. 
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7.3 The pseudo-real data setting 



Firstly we explain how the data have been created. We have considered for V the traded volumes 
of a very liquid security - namely the asset BNP - during an 11 day period. Then we selected the 
N most correlated assets (in terms of traded volumes) with the original asset. These assets are 
denoted Si, i = 1,...,N and we considered their traded volumes during the same 11 day period. 
Finally, the available volumes of each dark pool i have been modelled as follows using the mixing 
function 

/ EV ' 

VI < i < N, A := Pi 1(1- a t )V + a t S t 



'ESi 



where on, i 



1, . . . , N are the mixing coefficients, Pi, i = 1, . . . , N some scaling parameters and 
EV and ESi stand for the empirical mean of the data sets of V and Si. 



The shortage situation corresponds to A < 1 since it implies E £^ =] 
The simulations presented here have been made with four dark pools (N = 4). Since the data 
used here covers 11 days and it is clear that unlike the simulated data, these pseudo-real data are 
not stationary: in particular they are subject to daily changes of trend and volatility (at least). To 
highlight this resulting changes in the response of the algorithms, we have specified the days by 
drawing vertical doted lines. The dark pool pseudo-data parameters are set to 



Di 



< EV. 



P 



/0.1\ 

0.2 

0.3 
\0.2/ 



and 



a 



fOA\ 

0.6 

0.8 
\0.2/ 



and the dark pool trading (rebate) parameters are set to 

/0.01\ 
0.02 
9 ~ 0.04 
\0.06/ 

The mean and variance characteristics of the data sets of (V n ) n >i and (D") n >i, i = 1, ... ,4 
are the following: 





V 


D x 


D 2 


D 3 


D 4 


Mean 


955.42 


95.54 


191.08 


286.63 


191.08 


Variance 


2.01 x 10 e 


9.05 x 10 3 


4.29 x 10 4 


4.73 x 10 5 


5.95 x 10 4 



Firstly, we benchmarked both algorithms on the whole data set (11 days) as though it were 
stationary without any resetting (step, starting allocation, etc.). In particular, the running means 
of the performances are computed from the very beginning for the first 1500 data, and by a moving 
average on a window of 1500 data. As a second step, we proceed on a daily basis by resetting 
the parameters of both algorithms (the initial profit for the reinforcement algorithm (i.e. Ii = 0, 
1 < i < N) and the step parameter 7„ of the optimization procedure) at the beginning of every 
day. The performances of both algorithms are computed on each day. 



> Long-term optimization We observe that, except for the first and the fourth days where 
they behave similarly, the optimization algorithm is more performing than the reinforcement one. 
Its performance is approximately 30% higher on average (see Figure 3). 
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Relative Cost Reductions 




- Oracle 
Optimization 

- Reinforcement 



1 2 3 4 5 6 7 
Days 



9 10 11 



Performances 



- Oracle 
Optimization/Oracle 

- Reinforcement/Oracle 



1 2 3 4 5 6 7 
Days 



Figure 3: Long term optimization: Case N = 4, X^=i A < 1> < aj < 0.2 and = 1/N, 
l<i<N. 

This test confirms that the statistical features of the data are strongly varying from one day 
to another (see Figure [3]) , so there is no hope that our procedures converge in standard sense on a 
long term period. Consequently, it is necessary to switch to a short term monitoring by resetting 
the parameters of the algorithms on a daily basis as detailed below. 

> Daily resetting of the procedure We consider now that we reset each day all the parameters 
of the algorithm, namely we reset the step 7 n at the beginning of each day and the satisfaction 
parameters and we keep the allocation coefficients of the precedent day. We obtains the following 
results 



Relative Cost Reductions Performances 




1234567 89 10 11 1234567 8 9 10 11 

Days Days 



Figure 4: Daily resetting of the algorithms parameters: Case iV = 4, Y2i=i Pi < 1; < «j < 0.2 
and rf = 1/N 1 < i < N. 
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We observe (see Figure H|) that the optimization algoritm still significantly outperforms the 
reinforcement one, reaching more 95 % of the performance of the oracle. Furthermore, although 
not represented here, the allocation coefficients look more stable. 

8 Provisional remarks 

8.1 Toward more general mean execution functions 

One natural idea is to take into account that the rebate may depend on the quantity rV sent to 
be executed by the dark pool. The mean execution function of the dark pool can be modeled by 

Vr€[0,l], ip(r) = E(p(rV) mm(rV, D)) (8.1) 

where the rebate function p is a non-negative, bounded, non-decreasing right differentiable function. 
For the sake of simplicity, we assume that (V, D) satisfies (|3.4p . The right derivative of ip reads 

<p' r (r) = E (p' r (rV)Vmm (rV, D)) + E (p(rV)Vl {rV<D} ) , (8.2) 

with in particular <p'(0) = p(0) E(V1{£) >Q }) > as above. The main gap is to specify the function 
p so that ip remains concave which is the key assumption to apply the convergence theorem. 
Unfortunately the choice for p turns out to strongly depend on the (unknown) distribution of the 
random variable D. Let us briefly consider the case where V and D are independent and D has an 
exponential distribution £ (A). 

First note that the function g defined by g{u) := E(u AD), u > is given by 

1 _ e -«A 

Vu>0, g(u) = - 

so that, owing to the independence of V and D, 

Vr > 0, <p(r) = E{p{rV)g{rV)). 

At this stage, <p will be concave as soon as the function pg is so. Among all possible choices, 
elementary computations show that a class of possible choices is to consider p = g e with 0£ (0, A]. 
Of course this may appear as not very realistic since the rebate function is a structural feature of 
the different dark pools. 

However several numerical experiments not reproduced here testify that both algorithms are 
robust to a realistic choice for the function p e.g. a non-decreasing and stepwise constant. 

Another natural extension is to model the fact that the dark pool may take into account the 
volume rV to decide which quantity will really executed rather than simply the a priori deliverable 
quantity D. One reason for such a behaviour is that the dark pool may wish to preserve the 
possibility of future transactions with other clients. 

One way to model this phenomenon is to introduce a delivery function ip : IR^_ — > IR+, non- 
decreasing and concave w.r.t. its first variable and satisfying < ip(x,y) < y, so that the new 
mean execution function is as follows: 

<p(r) = pE(min(rV,ip(rV,D))) . (8.3) 
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It is clear that the function (p is concave (as the minimum of two concave functions) and bounded. 
In this case, the first (right) derivative of 92 reads 

ip' r (r) = pE(V (l{ rV< ^ rV)D)} +^' x (rV,D)l {rV >^ rViD)} )) (8.4) 

where ip' x denotes the right derivative with respect to x. In particular tp' r (Q) = pM(yijp >0 \) > 0. 

As concerns the implementations resulting from these new execution functions, the adaptation 
is straightforward. Note for the optimization procedure under constraints that, firstly, the "edge" 
functions Ri functions are not impacted by the type of the execution function. On the other hand 
the definition of the functions Hi or the updating of the variables I n for the reinforcement procedure 
should be adapted in accordance with the representations (I8.2D and (|8.4h of the new mean execution 
function ip. 

Example: We consider for modelling the quantity delivered by the dark pool i a function where 
we can define a minimal quantity required to begin to consum Dj, namely 

fpi(rV,Di) = A;l{rV> s ,D,} 
where Sj is a parameter of the dark pool i assumed to be deterministic. 



Pseudo-real data setting 



Relative Cost Reductions 




1 2 3 4 5 



9 10 11 



Figure 5: Shortage setting: Case N = 4, Yli=i ft < 1 5 < Oj < 0.2 and rf 
s = (0.3,0.2,0.2,0.3)*. 



1/N, 1 < i < N, 



8.2 Optimization vs reinforcement ? 

For practical implementation what conclusions can be drawn from our investigations on both proce- 
dures. Both reach quickly a stabilization/convergence phase close to optimality. The reinforcement 
algorithm leaves the simplex structurally stable which means the proposed dispatching at each time 
step is realistic whereas the stochastic Lagrangian algorithm in its present form sometimes needs 
to be corrected from time to time. This can be corrected by adding a projection on the simplex at 
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each step. We did not consider this variant from a theoretical point of view to keep our convergence 
proofs more elementary. 

In a high volatility context, the stochastic Lagrangian algorithm clearly prevails with perfor- 
mances that turn out to be significantly better. This optimization procedure also relies on estab- 
lished convergence results in a rather general framework (stationary a-mixing input data). However, 
given the computational cost of these procedures which is close to zero, a possible strategy is to 
implement them in parallel to get a synergistic effect. In particular, one may the reinforcement 
algorithm - which step parameter is structurally fixed equal to ^ - can be used to help tuning the 
constant c in the gain parameter 7„ = — of the stochastic Lagrangian. Doing so one may start with 
a small constant c, preventing a variance explosion of the procedure. Then based one may increase 
slowly this constant until the Lagragian outperforms the reinforcement procedure. 
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Appendix 

A Robbins-Zygmund with averaging innovation 

We consider an algorithm of the following form 

n +i = n - 7 „+i(G(0„, Y n ) + AM n+1 ), n > (A.l) 

where G is a Borel function from W d x W q to R d , {Y n ) n >o is a sequence of Revalued random vectors adapted 
to a filtration (J- n )n>i, &o is an J^-measurable M d -valued random vector independent of (Y n ) n >i, all defined 
on the same probability space J 7 , P), (AM„)„>i is a sequence of J-^-martingale increments and (7 n )n>i 
is a non-increasing sequence of positive real numbers going to as n goes to infinity. 

We will say that (Y n ) n >o is z/-averaging (under P) on a class of functions V + !P C P p {v) if, for every 
pe [l,+oo), 

V/ G V 0+ „, -V f(Y k ) — ► / fdv F-a.s. and in W(P). (A.2) 
n ^ Jr. 

Let f3€ (0, 1), let p e [1, oo). We denote by Vp tP the class of functions whose convergence rate in (|A.2[) P-a.s. 
and in L P (P) is 0(n- /3 ), namely 

V,, p = |/ : R"^ R | \ p f(Y k ) J fdv P ~ a - S - t 0(n -/>) J . (A.3) 

Now we are in a position to state the convergence theorem. 

Theorem A.l (A Robbins-Zygmund like Lemma) Let G : M. d x — > R d be a Borel function, let (Y n ) n >o 
be an T n -adapted v-averaging sequence of P 9 -valued random vectors and let (AM„) B >i be a sequence of 
J- n -martingale increments. Assume that there exists a continuously differ entiable function L : W d — > P + 
satisfying 

VL is Lipschitz continuous and |VL| 2 < C (1 + L) (A. 4) 

such that the function G satisfies the following local weak mean-reverting assumption: 

V6»e R d , \/ye P", (VZ(0) | G(0,y) - G(9*,y)) > 0. (A.5) 
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Suppose there exists f3£ (0,1), pG [l,oo) such that 

G(e*,-)€V 0tP . (A.6) 
Moreover, assume that G satisfies the following linear growth assumption 

V0eR d ,VyeR«, \G{9,y)\ < <p(y)(l + L(9)f* (A.7) 

and 

E(|AM„ +1 | 2 |J„) < V 2 (Y n )(l + L{6 n )) (A.8) 

where the function ip satisfies sup„ ||y(Y^)||n v < +00 (convention jjj = +00). 

Let 7 = ("f n )n>i be a non-negative, non-increasing sequence of gain parameters satisfying 

^7„ = +oo, n 1 -?^ — >0, and ^ fc 1 "' 3 max (7^, |A 7fe+ i |) < +00. (A.9) 

n>l fc>l 

Assume that 80 is J-Q-adapted and L(6q) < +00, P-a.s. Then, the recursive procedure defined by fA.l)) 
satisfies 

L{0 n ) ^ ioo < +00 and ^ 7 n+1 (Vi(^„) | G(9 n , Y n ) - G(d*,Y n )) < +00 a.s. 

n>0 
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